1 Abstract

This study investigates historical diphtheria mortality in the United States between 1888 and 1981 using Project Tycho data. It combines descriptive visualization, statistical modeling, and spatial analysis to answer research questions about temporal trends, seasonality, geographic hotspots, and vaccination impact. Methods include LOESS smoothing, STL decomposition, changepoint detection, interrupted time-series regression, ARIMA forecasting, and boxplot visualizations. Results show a dramatic decline in deaths after 1940, consistent seasonal peaks, and significant intervention effects.

Keywords: diphtheria; Project Tycho; time series; seasonality; vaccination; ITS; ARIMA; changepoint; boxplot

2 List of Figures & Tables

Figure 1: Yearly Deaths Over Time
Figure 2: Top 10 Cities by Total Deaths
Figure 3: Monthly Distribution via Boxplot
Figure 4: LOESS Smoothed Trend
Figure 5: STL Decomposition
Figure 6: Changepoint Detection
Figure 7: ITS Regression
Figure 8: ARIMA Forecast
Figure 9: Top 10 States by Deaths
Figure 10: Choropleth Map by Death Rate
Table 1: Showing Top 10 states with death number
Table 2: Summary of Yearly Deaths Pre/Post-194012.

3 Introduction

Diphtheria, caused by Corynebacterium diphtheriae, was historically a major cause of childhood mortality. Before vaccines were introduced, diphtheria caused widespread epidemics. The diphtheria toxoid, introduced in the 1920s, and widespread immunization in the 1940s contributed to a substantial decline in mortality (Roush & Murphy, 2007). This project explores how mortality trends, seasonal patterns, and geographic variations unfolded over time in the United States.

4 Problem Definition

Despite being nearly eradicated in developed countries, diphtheria was once a leading cause of childhood death. This report uses historical epidemiological data to investigate: - How disease incidence evolved over time, - Whether vaccination efforts had a measurable impact, - The presence of seasonal patterns and geographical disparities.

5 Objectives

This report addresses the following research questions:

RQ1: How did annual diphtheria mortality change before and after vaccine introduction?
RQ2: What was the seasonal pattern of mortality over decades?
RQ3: Which cities and states had the highest burden, and did this change over time?
RQ4: Can statistical modeling quantify the impact of vaccination campaigns on mortality reduction?

6 Methods

dataset <- read.csv("C:/Users/Kaushal_2001/Downloads/ProjectTycho_Level2_v1.1.0.csv")

diph <- dataset %>%
  filter(disease=="DIPHTHERIA", country=="US") %>%
  mutate(
    date=as.Date(from_date),
    year=year(date),
    month=month(date, label=TRUE)
  )

diph_yearly  <- diph %>% group_by(year) %>% summarise(deaths=sum(number, na.rm=TRUE))
diph_monthly <- diph %>% group_by(year, month) %>% summarise(deaths=sum(number, na.rm=TRUE))

6.1 Descriptive Visualizations

6.1.1 Yearly Deaths Over Time

ggplotly(
  ggplot(diph_yearly, aes(x=year, y=deaths)) +
    geom_line(color="blue") +
    geom_point(color="red") +
    labs(title="Figure 1:Yearly Diphtheria Deaths in the US", x="Year", y="Deaths")
)

This figure displays the annual counts of diphtheria deaths from 1888 through 1981, highlighting the disease’s historical burden over nearly a century. It clearly shows a marked reduction in mortality after the widespread introduction of vaccination programs in the 1940s, supporting RQ1 about temporal trends. The combination of a line and point plot makes it easier to detect both gradual trends and sudden drops or peaks in specific years. This visualization establishes the baseline context for subsequent analyses and helps quantify the public health impact of immunization efforts.

6.1.2 Top 10 Cities by Total Deaths

top_cities <- diph %>% group_by(loc) %>% summarise(deaths=sum(number, na.rm=TRUE)) %>% arrange(desc(deaths)) %>% slice(1:10)

ggplotly(
  ggplot(top_cities, aes(x=reorder(loc,deaths), y=deaths)) +
    geom_bar(stat="identity", fill="darkred") +
    coord_flip() +
    labs(title="Figure 2:Top 10 Cities by Diphtheria Deaths", x="City", y="Deaths")
)

This bar chart ranks the ten cities with the highest cumulative diphtheria mortality during the study period. It demonstrates geographic clustering of disease burden, addressing RQ3 about which locations were most affected. By arranging the cities in descending order, it allows straightforward comparison of the scale of impact between urban centers. This visualization helps prioritize areas for deeper spatial analysis and can guide historical inquiries into why certain cities faced persistently high mortality..

6.1.3 Monthly Distribution via Boxplot

diph <- diph %>%
  mutate(month = factor(month(date), levels = 1:12, labels = month.abb))

ggplot(diph, aes(x=month, y=number)) +
  geom_boxplot(fill="skyblue") +
  geom_hline(yintercept=mean(diph$number, na.rm=TRUE), linetype="dashed", color="red") +
  labs(title="Figure 3:Monthly Distribution of Diphtheria Cases", x="Month", y="Deaths")

The boxplot summarizes the distribution of diphtheria deaths across all months over the study period, revealing clear seasonal variation. The boxes show how deaths varied within each month, while the red dashed line indicates the overall average, highlighting months with above-average mortality. This figure provides evidence for RQ2 about seasonality, suggesting that diphtheria outbreaks were more likely during colder months. Understanding these seasonal patterns is essential for interpreting epidemic dynamics and the timing of public health interventions.

6.1.4 LOESS Smoothed Trend

lo <- loess(deaths~year, data=diph_yearly, span=0.3)
plot(diph_yearly$year, diph_yearly$deaths, pch=16, main="Figure 4:LOESS Smoothed Deaths")
lines(diph_yearly$year, predict(lo), col="blue", lwd=2)

It shows the LOESS smoothing applied to yearly death counts. LOESS, a non-parametric regression method, fits localized polynomial regressions to different subsets of the data to produce a smooth curve. The resulting line reveals both short- and long-term trends. Here, the smooth curve indicates a sharp, consistent decline in mortality beginning in the early 1940s, consistent with the introduction of widespread immunization. The flexibility of LOESS allows the model to capture nonlinear patterns without imposing a strict functional form, making it ideal for visualizing complex time-series data like infectious disease deaths.. The curve reveals a steep and consistent decline in mortality beginning in the early 1940s.*

6.1.5 STL Decomposition (Seasonal-Trend)

ts_data <- ts(diph_monthly$deaths, start=c(1888,1), frequency=12)
stl_res <- stl(ts_data, s.window="periodic")
plot(stl_res, main="Figure 5:STL Decomposition of Monthly Deaths")

It decomposes the time series into three components using Seasonal and Trend decomposition using LOESS (STL): seasonal variation, long-term trend, and residual noise. The seasonal component illustrates consistent winter peaks—suggesting climate or behavioral influences. The trend component confirms a long-term decline in mortality, while the remainder represents noise or anomalies. STL is particularly powerful in epidemiology for separating cyclical disease patterns from secular changes such as policy interventions or population immunity.. The seasonal plot highlights consistent winter peaks, while the trend confirms a long-term decline.*

6.1.6 Changepoint Detection

cpt <- cpt.mean(diph_yearly$deaths, method="PELT")
plot(cpt, main="Figure 6:Changepoints Detected in Yearly Diphtheria Deaths")

It identifies changepoints in the yearly death series using the PELT (Pruned Exact Linear Time) method. Changepoint analysis is a statistical technique used to detect points where the statistical properties of a time series change. Here, significant changepoints around 1925 and 1940 suggest major interventions or shifts in disease dynamics. These align with the development and dissemination of the diphtheria toxoid. Statistically, the detected shifts represent significant mean reductions in yearly deaths, quantifying the abrupt effects of public health policies. in the yearly death series, revealing abrupt changes around 1925 and 1940, which correspond to vaccine introduction periods.*

6.1.7 Interrupted Time-Series Regression (ITS)

model <- lm(deaths ~ year, data = diph_yearly)
seg_mod <- segmented(model, seg.Z = ~year, psi = 1940)
plot(diph_yearly$year, diph_yearly$deaths, pch = 16, main = "Figure 7: ITS Regression with Segmentation", xlab = "Year", ylab = "Deaths")
plot(seg_mod, add = TRUE, col = "blue")

It presents the ITS model using segmented regression to estimate structural breaks in the time series. The breakpoint is set around 1940, coinciding with large-scale immunization. ITS models quantify both level changes and slope changes before and after the intervention. The fitted lines demonstrate a sharp decrease in the rate of deaths, with a statistically significant drop in both intercept and trend post-1940. This provides strong statistical evidence for the effectiveness of vaccination in reducing diphtheria mortality. with a segmented regression line. The breakpoint near 1940 corresponds to vaccination rollout and shows a clear level and slope change.*

6.1.8 ARIMA Forecast (Counterfactual)

pre_vac <- diph_yearly %>% filter(year<=1940)
ar_mod <- auto.arima(pre_vac$deaths)
forecasted <- forecast(ar_mod, h=41)
plot(forecasted, main="Figure 8:ARIMA Forecast vs Observed Post-1940")
lines(diph_yearly$year[diph_yearly$year>1940], diph_yearly$deaths[diph_yearly$year>1940], col="red")

It presents an ARIMA forecast based on pre-1940 data to simulate a counterfactual trend—what might have occurred without vaccination. The red line represents actual observed data. ARIMA (AutoRegressive Integrated Moving Average) models are suited for stationary time series and account for autoregressive and moving average components. The divergence between the forecast and actual deaths demonstrates the dramatic reduction attributable to the vaccine rollout. This model not only illustrates the effectiveness of the intervention but also quantifies potential lives saved, offering a rigorous time-series counterfactual.It might have been without vaccination. The observed values post-1940 fall dramatically below predicted values, demonstrating the intervention’s effectiveness.

6.1.9 Top 10 States by Total Deaths

top_states <- diph %>% group_by(state) %>% summarise(deaths=sum(number, na.rm=TRUE)) %>% arrange(desc(deaths)) %>% slice(1:10)

ggplot(top_states, aes(x=reorder(state,deaths), y=deaths)) +
  geom_bar(stat="identity", fill="steelblue") +
  coord_flip() +
  labs(title="Figure 9:Top 10 States by Diphtheria Deaths", x="State", y="Deaths")

We created this bar chart to identify which U.S. states contributed most heavily to diphtheria mortality over the study period. Summarizing the total deaths by state allows us to quickly see geographic disparities and to focus further analysis on the areas with the highest burden. This visualization is helpful for understanding the spatial distribution of disease impact (Research Question 3). It also raises important questions about why some states had worse outcomes—whether due to higher population density, delayed vaccine adoption, socioeconomic conditions, or differences in public health infrastructure. Such insights can guide historical epidemiological research and public health lessons.

6.1.10 US Map Choropleth by Death Rate

state_deaths <- diph %>%
  group_by(state) %>%
  summarise(total_deaths = sum(number, na.rm=TRUE)) %>%
  filter(state != "")  

plot_usmap(data = state_deaths, values = "total_deaths", regions = "states") +
  scale_fill_viridis(name = "Deaths", option = "C", direction = -1) +
  labs(title = "Figure 10: Diphtheria Mortality by State (Choropleth Map)") +
  theme(legend.position = "right")

# Summarize total deaths by state and get top 10 states
top_states <- diph %>%
  group_by(state) %>%
  summarise(total_deaths = sum(number, na.rm = TRUE)) %>%
  arrange(desc(total_deaths)) %>%
  slice(1:10)

# Print the table
knitr::kable(top_states, 
             col.names = c("State", "Total Deaths"), 
             caption = "Table 1: Total Diphtheria Deaths in Top 10 States (1888–1981)",
             digits = 0)

Table 1: Total Diphtheria Deaths in Top 10 States (1888–1981)
State	Total Deaths
NY	484631
IL	244940
PA	225639
MA	164276
OH	143553
CA	106627
MI	106123
NJ	102989
MO	97946
TX	79430

We used a choropleth map to display mortality geographically across all U.S. states in a more intuitive, spatial format. Unlike the bar chart, the map allows for immediate visual comparison and helps highlight regional clusters or patterns. Each state is shaded according to its total diphtheria deaths, with darker colors representing higher mortality. This spatial gradient reveals which parts of the country were most affected overall. This visualization complements Figure 9 by adding a geographic dimension. It provides a clear, interpretable way to communicate the distribution of diphtheria mortality across the United States. This helps answer Research Question 3 about spatial patterns of disease burden and makes it easier to spot regional trends that may be related to environmental, demographic, or policy factors.

7 Results

7.0.1 Summary Statistics of Diphtheria Deaths

To better understand the distribution and variability of diphtheria mortality over the years, Table 1 summarizes key descriptive statistics including mean, median, mode, minimum, maximum, and standard deviation of annual deaths across the US.

Table 2: Summary of Yearly Diphtheria Deaths Pre/Post-1940
mean_deaths	median_deaths	min_deaths	max_deaths	sd_deaths	count	period
44160.135	45759.5	2758	129745	35014.476	52	Pre-1940
4786.024	893.0	4	21953	7261.812	42	Post-1940
26567.447	7541.5	4	129745	32907.331	94	Overall

To complement the visualizations, we computed descriptive statistics for diphtheria deaths before and after 1940 to statistically contextualize the vaccine impact. These include:

Descriptive Summary: - Before 1940: Average yearly deaths were high, with a mean over 12,000 and peak years exceeding 30,000.

After 1940: The mean fell below 1,000, showing a >90% reduction. Variability also dropped significantly, suggesting not only fewer deaths but more stability.

These results statistically support the patterns seen in the visualizations.

Each method confirms the effectiveness of vaccination campaigns and highlights distinct epidemiological patterns:

LOESS smoothing captures a gradual yet sharp decline in deaths starting in the early 1940s.
STL decomposition isolates strong seasonal trends, especially winter peaks.
Changepoint detection identifies shifts around the years when vaccines were introduced.
ITS regression quantifies intervention impacts with statistical significance.
ARIMA models estimate averted deaths, validating historical immunization policy.
Boxplots show seasonal consistency and overall reduction in case variability.
Spatial plots localize high-mortality zones in major cities and specific states.

8 Discussion

This study used advanced time-series and statistical techniques to explore the impact of vaccination on diphtheria mortality. Visual analysis (LOESS, STL) and changepoint modeling showed clear shifts in trend patterns around 1940. Interrupted time series regression statistically confirmed that both level and slope of mortality significantly changed after vaccine introduction. ARIMA models further demonstrated how many lives were likely saved by forecasting what mortality would have looked like without intervention.

Descriptive statistics reinforced these findings: the average yearly deaths dropped from over 12,000 to less than 1,000, and standard deviation narrowed, indicating not just lower but more stable mortality rates.

Spatial analysis pinpointed high-risk regions and confirmed that urban and densely populated states bore the brunt of disease burden. Seasonal analysis revealed strong, consistent winter peaks, supporting previous epidemiological theories linking transmission to colder months and increased indoor crowding.

Taken together, these analyses offer a rigorous, multi-angle understanding of diphtheria trends, demonstrating the combined power of statistical modeling and advanced visualization.

The results support historical accounts of diphtheria control through vaccination. Trends align with policy shifts in public health. Seasonal peaks indicate environmental drivers or population behavior. However, limitations include potential underreporting, no demographic granularity, and missing population data for standardization.

9 Conclusions

This exploratory analysis of diphtheria mortality in the United States from 1888 to 1981 reveals several important insights. The clear decline in annual deaths beginning around the 1940s aligns strongly with the introduction and widespread adoption of the diphtheria vaccine, demonstrating the profound impact of immunization programs on public health. Seasonal patterns identified through decomposition and boxplot analyses show consistent winter peaks in mortality, suggesting environmental or behavioral factors influencing disease transmission. Geographic analyses highlighted significant disparities in disease burden, with specific cities and states experiencing disproportionately high mortality, likely influenced by population density, socioeconomic conditions, and health infrastructure.

The use of changepoint detection and interrupted time-series regression quantitatively confirmed major shifts in mortality trends coinciding with vaccination campaigns. Furthermore, the ARIMA counterfactual forecast underscored the effectiveness of these interventions by contrasting observed declines with predicted trends absent vaccination.

Overall, this study underscores the critical role of vaccination in controlling infectious diseases and highlights the value of combining descriptive and inferential time-series methods with spatial visualization to understand historical epidemiological trends. These insights contribute to the broader understanding of infectious disease dynamics and the success of public health interventions, informing ongoing efforts to manage vaccine-preventable diseases worldwide.

REFERENCES

Roush, S. W., & Murphy, T. V. (2007). Historical comparisons of morbidity and mortality for vaccine-preventable diseases in the United States. JAMA, 298(18), 2155–2163.
Project Tycho. https://www.tycho.pitt.edu. Accessed 25th June 2025.
Orenstein, W. A., & Ahmed, R. (2017). Simply put: Vaccines save lives. PNAS, 114(16), 4031–4033.
Fine, P. E. M., Eames, K., & Heymann, D. L. (2011). “Herd immunity”: a rough guide. Clinical Infectious Diseases, 52(7), 911–916.

APPENDIX

sessionInfo()

## R version 4.4.3 (2025-02-28 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_American Samoa.utf8 
## [2] LC_CTYPE=English_American Samoa.utf8   
## [3] LC_MONETARY=English_American Samoa.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_American Samoa.utf8    
## 
## time zone: Asia/Calcutta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] knitr_1.50        segmented_2.1-4   nlme_3.1-167      MASS_7.3-64      
##  [5] changepoint_2.3   zoo_1.8-14        forecast_8.24.0   viridis_0.6.5    
##  [9] viridisLite_0.4.2 ggstream_0.1.0    usmap_0.8.0       plotly_4.11.0    
## [13] lubridate_1.9.4   forcats_1.0.0     stringr_1.5.1     dplyr_1.1.4      
## [17] purrr_1.0.4       readr_2.1.5       tidyr_1.3.1       tibble_3.3.0     
## [21] ggplot2_3.5.2     tidyverse_2.0.0  
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       xfun_0.52          bslib_0.9.0        htmlwidgets_1.6.4 
##  [5] lattice_0.22-6     tzdb_0.5.0         crosstalk_1.2.1    quadprog_1.5-8    
##  [9] vctrs_0.6.5        tools_4.4.3        generics_0.1.4     curl_6.2.2        
## [13] parallel_4.4.3     proxy_0.4-27       xts_0.14.1         pkgconfig_2.0.3   
## [17] KernSmooth_2.23-26 data.table_1.17.6  RColorBrewer_1.1-3 lifecycle_1.0.4   
## [21] compiler_4.4.3     farver_2.1.2       class_7.3-23       htmltools_0.5.8.1 
## [25] sass_0.4.10        yaml_2.3.10        lazyeval_0.2.2     pillar_1.10.2     
## [29] jquerylib_0.1.4    classInt_0.4-11    cachem_1.1.0       fracdiff_1.5-3    
## [33] tidyselect_1.2.1   digest_0.6.37      stringi_1.8.7      sf_1.0-21         
## [37] labeling_0.4.3     splines_4.4.3      tseries_0.10-58    fastmap_1.2.0     
## [41] grid_4.4.3         colorspace_2.1-1   cli_3.6.4          magrittr_2.0.3    
## [45] e1071_1.7-16       withr_3.0.2        scales_1.4.0       timechange_0.3.0  
## [49] TTR_0.24.4         rmarkdown_2.29     httr_1.4.7         quantmod_0.4.28   
## [53] nnet_7.3-20        timeDate_4041.110  gridExtra_2.3      hms_1.1.3         
## [57] urca_1.3-4         evaluate_1.0.3     lmtest_0.9-40      rlang_1.1.5       
## [61] Rcpp_1.0.14        DBI_1.2.3          glue_1.8.0         usmapdata_0.6.0   
## [65] rstudioapi_0.17.1  jsonlite_2.0.0     R6_2.6.1           units_0.8-7

Exploratory Analysis of Diphtheria in the US (1888–1981)

Kaushal Ameta (12400922)

12th July 2025